function [r] = dirichlet_sample(a,m,n)

% DIRICHLET_SAMPLE   Sample from Dirichlet distribution.
% DIRICHLET_SAMPLE(a) returns a probability vector sampled from a 
% Dirichlet distribution with parameter vector A.
% DIRICHLET_SAMPLE(a,n) returns N samples, collected into a matrix, each 
% vector having the same orientation as A.
% Copyright: Ayan Acharya

if nargin < 3
  n = 1;
end

row = (size(a,1) == 1);

a = a(:)*m;
y = gamrnd(repmat(a, 1, n),1);
r = sum(y,1);
r(find(r == 0)) = 1;
r = y./repmat(r, size(y,1), 1);

% % %% verifying frequency
% % 
% % [~,rmaxI] = max(r);
% % hist(rmaxI,unique(rmaxI));

if row
  r = r';
end

end

